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Abstract. Correlated fermions are of high interest in condensed matter (Fermi 
liquids, Wigner molecules), cold atomic gases and dense plasmas. Here we propose a 
novel approach to path integral Monte Carlo (PIMC) simulations of strongly degenerate 
non-ideal fermions at finite temperature by combining a fourth-order factorization of 
the density matrix with antisymmetric propagators, i.e., determinants, between all 
imaginary time slices. To efficiently run through the modified configuration space, we 
introduce a modification of the widely used continuous space worm algorithm, which 
allows for an efficient sampling at arbitrary system parameters. We demonstrate 
how the application of determinants achieves an effective blocking of permutations 
with opposite signs, leading to a significant relieve of the fermion sign problem. 
To benchmark the capability of our method regarding the simulation of degenerate 
fermions, we consider multiple electrons in a quantum dot and compare our results 
with other ab initio techniques, where they are available. The present permutation 
blocking path integral Monte Carlo approach allows us to obtain accurate results even 
for = 20 electrons at low temperature and arbitrary coupling, where no other ab 
initio results have been reported, so far. 
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1. Introduction 

The ab initio simulation of strongly degenerate nonideal fermions at finite temperature 
is of high current importance for many fields. The numerous physical applications 
include electrons in a quantum dot [1, 2, 3, 4], fermionic bilayer systems [5, 6, 7], the 
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homogeneous electron gas [8, 9, 10], dense two-component plasmas [11, 12, 13] in stellar 
interiors and modern laser compression experiments (warm dense matter) [14, 15] and 
inertial fusion [16]. Despite remarkable recent progress, existing simulation methods 
face serious problems. 

The widely used path integral Monte Carlo (PIMC) method, e.g. [17], is a 
highly successful tool for the ab initio simulation of both distinguishable particles 
(“boltzmannons“, e.g. [18, 19]) and bosons [17] and allows for the calculation of quasi- 
exact results for up to ~ 10^ particles [20] at hnite temperature. However, the 
application of PIMC to fermions is hampered by the notorious sign problem [21], which 
renders even small systems unfeasible for state of the art techniques and has been 
revealed to be iVP-complete for a given representation [22]. With increasing exchange 
effects, permutation cycles with opposite signs appear with nearly equal frequency and 
the statistical error increases exponentially. For this reason, standard PIMC is applicable 
to fermions only at weak degeneracy, that is, at relatively high temperature or low 
density. 

The recently introduced conhguration path integral Monte Carlo (CPIMC) method 
[23, 24, 9] exhibits a complementary behavior. This conceptually different approach can 
be interpreted as a Monte Carlo simulation on a perturbation expansion around the 
ideal quantum system and, therefore, CPIMC excells at weak nonideality and strong 
degeneracy. Unfortunately, the physically most interesting region, where both fermionic 
exchange and interactions are strong simultaneously, remains out of reach. 

A popular approach to extend standard PIMC to higher degeneracy is Restricted 
PIMC (RPIMC) [25], also known as hxed node approximation. This idea requires 
explicit knowledge of the nodal surfaces of the density matrix, which are, in general, 
unknown and one has to rely on approximations, thereby introducing an uncontrollable 
systematic error. In addition, it has been shown analytically [26, 27] that RPIMC does 
not reproduce the exact density matrix in the limit of the ideal Fermi gas and, therefore, 
the results become unreliable at increasing degeneracy [9]. 

Recently, DuBois et al. [28] have suggested that, at least for homogeneous systems, 
the individual exchange probabilities in PIMC are independent of the conhguration of 
other permutations present and that permutation frequencies of large exchange cycles 
can be extrapolated from few-particle permutations. This would allow for a signihcant 
reduction of the conhguration space and a drastic reduction of the sign problem. While 
hrst simulation results with this approximation for the short-range interacting ^He 
are in good agreement with experimental data [28], the existing comparison [9] for 
long-range Coulomb interaction is insufficient to assess the accuracy and, in addition, 
inhomogeneous systems remain out of reach. 

Another possibility to relieve the sign problem in fermionic PIMC without 
introducing any approximations is the usage of antisymmetric imaginary time 
propagators, i.e., determinants [29, 30, 10, 31]. It is well known that the sign problem 
becomes more severe with an increasing number of propagators arising from the Trotter- 
type factorization of the density operator. Consequently, it has been proposed to 
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Figure 1. Illustration of the capability of PB-PIMC - In panel (A), the average sign 
S from different methods is plotted versus the coupling parameter A, Eq. (31), for 
N = 20 electrons in a quantum dot at /3 = 3.0 (oscillator units). Region [I] denotes the 
weakly nonideal Fermi gas, [II] the transition region and [III] the strongly correlated 
regime. CPIMC (PIMC) is limited to weak (strong) coupling, i.e. to the region left 
(right) of the blue (green) line. Panel (B) shows a comparison of density profiles n(r), 
plotted versus the distance to the center of the trap r, across the entire coupling range. 


combine the antisymmetric propagators with a higher order factorization [32, 33, 34, 35] 
of the density matrix. This has recently allowed to obtain an accurate estimate of the 
ground state energy of degenerate, strongly nonideal electrons in a quantum dot [36]. 

In the present work, we extend this idea to finite temperature. For this purpose, we 
combine a fourth-order propagator derived in [37], which has already been succesfully 
applied to PIMC by Sakkos et al. [38], with a full antisymmetrization on all time slices 
to simulate fermions in the canonical ensemble. We demonstrate that the introduction 
of determinants effectively allows for the combination of N\ configurations from usual 
PIMC into a single configuration weight, thereby reducing the complexity of the problem 
and blocking both positive and negative weights to drastically increase the sign. To 
efficiently exploit the resulting configuration space with the Metropolis algorithm [39] 
at arbitrary parameters, we develop a set of Monte Carlo updates similar to the usual 
continuous space worm algorithm (WA) [20, 40]. 

To demonstrate the capability of our permutation blocking path integral Monte 
Carlo (PB-PIMC) method, we consider Coulomb interacting fermions in a 2D harmonic 
conhnement, cf. Eq. (30), which can be experimentally realized e.g. by spin-polarized 
electrons in a quantum dot [1, 2, 3, 4]. Figure 1 (A) shows the average sign S' for 77 = 20 
electrons, plotted versus the coupling strength A, cf. Eq. (31). CPIMC is applicable in 
the weakly nonideal regime [I], where the system is predominantly shaped by the Fermi 
statistics. In contrast, standard PIMC allows one to accurately simulate systems in the 
strongly coupled regime [III], where exchange effects are not yet dominating, and bosons 
and fermions exhibit a very similar behavior. The PB-PIMC method, as will be shown 
in this work, is applicable over the entire coupling range yielding reasonably accurate 
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results with acceptable computational effort. Interestingly, this includes the physically 
most interesting transition region [11], where both the Coulomb repulsion and quantum 
statistics govern the system. Here no ab initio results have been reported to this date, 
except for very small particle numbers, since PIMC and CPIMC fail, due to the sign 
problem. In panel (B), we show density profiles from all three regimes. Evidently, the 
transition from the strongly coupled system with a pronounced shell structure (A = 15) 
to the nearly ideal Fermi gas with the characteristic weak density modulations (A = 0.1) 
can be resolved. 

In the remainder of this work, we introduce the PB-PIMC method in detail. We 
show that the optimal choice of two free parameters of the fourth-order factorization 
allows for a calculation of energies and densities with an accuracy of the order of 0.1% 
with as few as two or three propagators, even in the low temperature regime. We 
calculate energies and densities from PB-PlMC for iV = 20 electrons at low temperature 
over the entire coupling range. We find excellent agreement with both PIMC and 
CPIMC in the limitting cases of strong and weak coupling, respectively, and perform 
simulations in the transition regime, where no other ab initio results are available. 
Finally, we investigate the performance behavior of our method when the system size is 
varied. 


2. Theory 


2.1. Idea of permutation blocking path integral Monte Carlo 


We consider the canonical ensemble (the particle number N, volume V and inverse 
temperature /3 = l/ksT are fixed) and write the partition function in coordinate 
representation as 

^ ^ ’ ( 1 ) 
ctSSjv 

where R = {ri,...,rAr} contains the coordinates of all particles and denotes the 
exchange operator corresponding to a particular element a from the permutation group 
S'tv. The Hamiltonian is given by the sum of the kinetic (K) and potential (V) energy, 
H = K + V. For the next step, we use the group property of the density operator 

P = = , ( 2 ) 

with e = fl/P, and insert P — 1 unities of the form 1 = / dR^ |Ra) (Rq,|. This gives 
Z = f dRo ... dRp_i U ( ®Sn(cr) (R„| |7r^Ro+i) j . (3) 

a=0 \ ■ aeSN J 


Note that we have exploited the permutation operator’s idempotency property in Eq. 
(3) to introduce antisymmetry on all P imaginary time slices. Following Sakkos et al. 
[38], we introduce the factorization from [37], 


e 


-eH 


g-tiek'g-V2eW'i_2ai ^-HeK^-v^eWa^ ^-2toeK 




(4) 
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for each of the exponential functions in Eq. (3). By including double commutator terms 
of the form 

N 






( 8 ) 


2=1 


we have to evaluate the total force on each particle, Fj = —ViF(R), and Eq. (4) is 
accurate to fourth order in e. The explicit form of the modihed potential terms W is 
given by 


Wa, =V + -aie^l-^jF 
Vi \ m ^^ 




N 


and 


Vl 


= E + ^(1 - ai)e2 I - |F. 

Vo \ m 


( 6 ) 


2=1 




N 


V2 


2=1 


There are two free parameters in Eq. (4), namely 0 < oi < 1, which controls the relative 
weight of the forces on a particular slice, and 0 < to < (1 ~ l/\/3)/2, which determines 
the ratio of the, in general, non-equidistant time steps between ” daughter “ slices, cf. 
Fig. 2. All other factors are calculated from these choices: 


“0 = yx 1 - 


12 


+ 


1 - 2to 6(1 - 2tof 


Vl = 


6(1 - 2to)2 ’ 

^2 = 1 — 2ni and 
_ 1 

— X ~ ^0 


( 7 ) 


The fourth-order approximation of the imaginary time propagator is visualized in 
Fig. 2. The inverse temperature (3 has been split into P = 4 intervals of length e, which 
are further divided into three, in general, non-equidistant sub-intervals. Thus, for each 
main ”bead“ r^, there exist two daughter beads, and TaB- 

Let us for a moment ignore the antisymmetry in Eq. (3) and evaluate the imaginary 
time propagator in a straightforward way [38]: 


(Role "-^iRa+i) 


dRo-AdRo-s 





a 


N 

Pa(h *)Pq:,a(L *)Pas(L ) 

i=l 

with the dehnitions of the potential terms 


( 8 ) 


Va = niE(Ra) -f voV (Rqa) + viV (R^b) , 

N 

Fa= (Qi|Fq,^;|^ -|- (1 — 2ai) |FQ,A,i|^ + Ql|FaB,^|^) 5 

2=1 


(9) 
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Figure 2. Illustration of the configuration space - In the left panel, the imaginary 
time is plotted versus the (arbitrary) spatial coordiante x. Each time step of length e 
is further divided into three non-equidistant subintervals, with two ’’daughter" slices 
A and B. The right panel illustrates the combination of all 3PNI possible trajectories 
into a single configuration weight IF(X). Between each two adjacent time slices, both 
the connection between beads from the same particle (diagonal elements of the diffusion 
matrix, the blue and red lines) and between beads from different particles (off-diagonal 
elements, the green lines) are efficiently grouped together to improve the average sign. 


and the diffusion matrices 

Pa{i,j) =A“fexp 


— (r ■ 

tie 


X 


W,i)' 


PaA{i,j) = At^fexp , 

PaB{i,j) = A^ijexp , 


( 10 ) 


where Xg denotes the thermal wavelength A^ = 2nh‘^/3/m and D is the dimensionality of 
the system. Thus, the matrix elements of Eq. (10) are equal to the free particle density 
matrix, Pa{i,j) = Po{pa,ji'^aA,iiic)- The permutation operator commutes with both p 
and H and we are, therefore, allowed to artihcially introduce the antisymmetrization 
between all 3P slices without changing the result. This transforms Eq. (8) to 

1 


N\ 


^ sgn((T) (Ro 


^ |7r.R„+i) = ( 


a-£Sj^ 


N\) 


d A d-R< q: 


e " ““-■^“det(p„)det(p„yi)det(p„s) 


Finally, this gives the partition function 

. P-i 


Z = 


(TV!) 


3P 


dxfl 




det (p«) det (p„A) det (p„B) 


(11) 


( 12 ) 


Q = 0 


and the integration is carried out over all coordinates on all 3P slices: 


dX — dRg ... dRp_idRoyi... dRp—iy^dRop ... dRp_ip 


(13) 
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The benefits of the partition fnnction Eq. (12) are illnstrated in the right panel 
of Fig. 2 where the beads of two particles are plotted in the r-x-plane. In the usnal 
PIMC formnlation (withont the determinants), each of the particles would correspond 
to a single closed trajectory as visualized by the blue and red connections. To 
take into account the antisymmetry of fermions, one would also need to sample all 
conhgurations with the same positions of the individual beads but different connections 
between adjacent time slices, which have both positive and negative weights. By 
indroducing determinants between all slices, we include all N\ possible connections 
between beads on adjacent slices (the green lines) into a single conhguration weight 
and the usual interpretation of mapping a quantum system onto an ensemble of 
interacting ringpolymers [41] is no longer appropriate. Therefore, a large number of sign 
changes, due to different permutations, are grouped together resulting in an efficient 
compensation of many terms (blocking), and the average sign (cf. Eq. (22)) in our 
simulations is signihcantly increased [31]. 


2.2. Energy estimator 


The total energy E follows from the partition function via the familiar relation 

—■ 

Substituting the expression from Eq. (12) into (14) and performing a lengthy but 
straightforward calculation gives the hnal result for the thermodynamic (TD) estimator 

P-l N N 

\ePXl, 

k = 0 K=1 ^=1 \ 


E = 


3DN 


2e 




^kA,^y 


- ^{rkA,> 


ePX 
1 


tie 

P-l 


P 


k=0 








{rkB,K - rk+uf 


(15) 


( 14 + 3e^uo—Fk 


with the dehnitions 

= {pkl).^iPkA)^, 

To split the total energy into a kinetic and a potential part, we evaluate 

_m 

PZ dm 


( 16 ) 


(17) 
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and find the TD estimator of the kinetic energy 

P-l N N 


fe = 0 K=l g=l . 


2e 


-- ^kA,{) + ^py2 {^kA,K - 


ePX 


tie 


tie 


TTVl/^f 


(rfcB,Ac - rfc+i,^)' 


p-i 




k=0 


k? ~ 
e Uo—Fk 
m 


Thns, the estimator of the potential energy is given by 


V = E-K = -y^(Vu + 2Cuo-h 

P \ m 

k=0 


(18) 


(19) 


We notice that the forces contribnte to both the kinetic and the potential energy. For 
completeness, we mention that, for an increasing nnmber of propagators, P —)• cx), 
the first and second terms in Eq. (15) diverge, which leads to a growing variance and, 
therefore, statistical nncertainty of both E and K. To avoid this problem, one might 
derive a virial estimator, e.g. [42], which reqnires the evaluation of the derivative of 
the potential terms instead. However, since we are explicitly interested in performing 
simulations with few propagators to relieve the fermion sign problem, the estimator 
from Eq. (15) is sufficient. 


3. Monte Carlo algorithm 


In section 2, we have derived an expression for the partition function Z, Eq. (12), which 
incorporates determinants of the diffusion matrices between all 3P time slices, thereby 
combining 3PN\ different configurations from the usual PIMC into a single weight 
fF(X). However, each determinant can still be either positive or negative, depending 
on the relative magnitude of diagonal and off-diagonal elements. Hence, we apply the 
Metropolis algorithm [39] to the modified partition function 

Z' = ydX|W(X)| , (20) 


and calculate fermionic expectation values as 

(OS)’ 


{ 0 ),= 


(S) 


with the definition of the average sign 

(S)' = /dx |fr(x)|s(x) 


and the signum of the conhguration X, 


( 21 ) 


( 22 ) 


P-l 

*S'(X) = [sgn(det(p„))sgn(det(p«^))sgn(det(paB))] • (23) 

Let us summarize some important facts about the configuration space defined by Eq. 

( 20 ): 
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Figure 3. Influence of the imaginary time step e on the efficiency of the permutation 
blocking - Two configurations of = 2 particles are visualized in the r-a:-plane. In the 
left and right panel, there are P = 2 and P = 5 time slices, respectively (daughter slices 
are neglected for simplicity). Only with few propagators, the thermal wavelength Ag of 
a single propagator is comparable to the mean interparticle distance d, which is crucial 
for an efhcient grouping of permutations into a single configuration weight. With 
increasing P, diagonal (red and blue lines) and off-diagonal (green lines) distances are 
no longer of the same order and the permutation blocking is inefficient. 


(i) With increasing number of propagators P, the effect of the blocking decreases and, 
for P —>■ oo, the sign converges to the sign of standard PIMC. Blocking is maximal 
if and X 2 toe are comparable to the average interparticle distance d, cf. Fig. 3. 
Only in such a case, there can be both large diagonal and off-diagonal elements in 
the diffusion matrices. 

(ii) Conhguration weights |fF(X)| can only be large, when at least one element in each 
row of each diffusion matrix is large. Therefore, we sample either large diagonal or 
large off-diagonal elements. Blocking happens naturally as a by-product and does 
not have to be specihcally included into the sampling. This also means that we 
have to implement a mechanism to sample exchange, i.e., to switch between large 
diagonal and off-diagonal diffusion matrix elements. 

(hi) There are no hxed trajectories. Therefore, beads do not have a previous or a 
next bead, as in standard PIMC. For an efhcient and flexible sampling algorithm, 
we temporarily construct artihcial trajectories and choose the included beads 
randomly. 

The most efhcient mechanism for the sampling of exchange cycles in standard PIMC 
is the so-called worm algorithm [20, 40], where macroscopic trajectories are naturally 
realized by a small set of local updates which enjoy a high acceptance probability. In the 
rest of the section, we modify this algorithm to be applicable to the new conhguration 
space without any hxed connections between individual beads. 
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Figure 4. Illustration of the sampling scheme (left) and the extended configuration 
space (right) - In the left panel, an artificial trajectory (pink curve) with four missing 
beads is plotted in the r-x-plane. The new coordinates (green circles) are sampled 
according to a Gaussian (blue curves) around the intersection of the connecting straight 
lines between the previous and last bead with the current time slice (black crosses). The 
right panel gives an example for an open configuration in the extended configuration 
space with two special beads which are denoted as ”head“ and ”tail“. There are 
only N — 1 beads on eight time slices, going forward in imaginary time starting from 
Thead = r 2 A- The circles, triangles and squares distinguish beads from three different 
particles and the empty symbols at the right boundary indicate the missing beads on 
a particular slice. 


3.1. Sampling scheme 


To take advantage of the main benefits from the usual continuous space worm algorithm, 
we will temporarily construct artihcial trajectories and sample new beads according to 
standard PIMC techniques, e.g. [43]. The initial situation for our considerations is 
illustrated in the left panel of Fig. 4, where a pre-existing trajectory (pink curve) with 
four missing beads in the middle is shown in the r-x-plane. We choose the sampling 
probability to close the conhguration as 




sample 


11^=0 7 I ’ z + 1 : ^^ 2+1 ' 7~0 


(24) 


Po(ro5 I’m, 'Tm — 'T’o) 

which results in the consecutive generation of M — 1 new coordinates r^, i G [1, M — 1], 
according to 






PoiPi-1) ^M,Tm — P-l) 



V 2a2 ) 


-P) 


1 


(25) 


which is a Gaussian (cf. the blue curves in Fig. 4) with the variance 

2 {n - Ti-i){TM - Ti) 

^ 0 •! 


m 


TM — Ti-1 


(26) 
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around the intersection of the connection between the previous coordinate, ri_i, with 
the end point ym and the time slice r* 


'^i—1 

-fi-i H-rju 

TM — Ti-i TM — Ti-i 


(27) 


3.2. Artificial worm algorithm 


In the usual WA-PIMC, the configuration space is defined by the Matsubara Green 
function (MGF, e.g. [44]) which implies that the algorithm does not only allow for the 
change of the particle number N (grand canonical ensemble) but, in addition, requires 
the generation of configurations with a single open path, the so-called worm. However, 
in the PB-PIMG configuration space defined by Eq. (12), there are no trajectories and, 
therefore, no direct realization of a worm is possible. Instead, we consider an extended 
ensemble, which combines closed configurations with a total of 3iVP beads and open 
configurations, where on some consecutive time slices the number of beads is reduced 
by one, to iV — 1. Such a configuration is illustrated in the right panel of Fig. 4. 
There are two special beads which are denoted as ”head“ and ”taiP‘ and the triangles, 
circles and squares symbolize beads from three different particles. There are eight beads 
from different particles missing (indicated by the empty symbols at the right boundary) 
between Thead = 'T 2 A and rtaii = ti^, going forward in imaginary time. 

For most slices, the computation of the diffusion matrix allows for no degree of 
freedom in the extended ensemble. We define the latter in a way, that the head bead 
does not serve as a starting point for the elements but is treated as if it was missing. This 
is justihed because, otherwise, there does not necessarily exist a large matrix element in 
this particular row because no artificial connection has been sampled on the next slice. 
For the configuration from Fig. 4, the diffusion matrix of the head’s time slice is given 
by 


P2A 


det(p: 




'^Po(ri,2A,ri,2s,tie) Po(ri,2A, r2,2B, He) 0^ 
= 1 11 
\Po(r3,2A,ri,2B,He) Po(r3,2A, r2,2B, tie) Oj 

= det I ’’1-25) He) Po(ri, 2 A, ^ 2 , 23 , he) 

\Po(r3,2A, ri, 2 B, tie) Po(r3,2A, ^ 2 , 23 , he) 


(28) 


All diffusion matrices with At—1 beads on their slices are computed in the same way. The 
other degree of freedom for which the extended ensemble allows is the choice whether 
the tail will be included as the final coordinate in the diffusion matrix or not. Here, it 
makes sense to allow for this possibility, because there does exist at least a single large 
element in this particular row anyway. The corresponding matrix for the configuration 
from Fig. 4 looks like 


P2A — 


'^Po(ri,i,ri,iA,He) Po(ri,i,r 2 ,iA,He) Po(ri,i, r3,iA, He)^ 

Po(r2,i,ri,iA,He) Po(r2,i,r2,iA,He) Po(r2,i, r3,iA, He) 

VI 1 1 / 


(29) 
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Table 1. Convergence of the energy for TV = 4, A = 1.3 and j3 = 5.0 for selected 
parameter combinations shown in Fig. 5. 


Simulation 

E 

V 

K 

S 

^P = 2 

12.1870(8) 

9.0157(8) 

3.1713(9) 

0.4950(5) 

'^P = 2 

12.3188(7) 

9.0840(6) 

3.2348(6) 

0.3790(5) 

^P= 15 

12.292(7) 

9.084(1) 

3.207(7) 

0.02456(8) 

15 

12.294(4) 

9.0827(9) 

3.214(4) 

0.01911(4) 

CPIMC 

12.293(3) 

- 

- 

- 


^ to = 0.04, ai = 0.0 
to = 0.13, ai = 0.33 


However, we emphasize that the particular choice of the extended ensemble does not 
influence the extracted canonical expectation values as long as detailed balance is 
fulhlled in all updates. We have developed a simulation scheme which consists of four 
different types of moves that ensure detailed balance and ergodicity. The updates are 
presented in detail in Appendix A. 


4. Simulation results 


As a test system to benchmark our method, we consider N spin-polarized electrons 
in a quantum dot [1, 2, 3, 4], which can be described approximately by a harmonic 
conhnement with a frequency hi. We use oscillator units, i.e., the characteristic energy 
scale Eq = hVt and oscillator length I = ^Jh/Vtm, and obtain the dimensionless 
Hamiltonian 


N 




2 ^ 

i=l 


N 


2 = 1 


N 

E 

Kj 


A 


Yi - r. 


(30) 


with the coupling parameter 
^ ^ 

being defined as the ratio of Coulomb and oscillator energy. For large A, the electrons are 
strongly coupled and exchange effects become negligible (region [HI] in Fig. 1), while, 
for A -C 1, the ideal Fermi gas will be approached and the system is governed by the 
fermionic exchange (region [I] in Fig. 1). To confirm the quality of our simulations, 
we compare the results at weak and strong coupling with CPIMC and standard PIMC, 
respectively, where they are available. 


4 -1. Optimal choice of ai and to 

We start the discussion of the simulation results by investigating the effects of the two 
free parameters oi and to on the convergence of two different observables, namely the 
energy E and radial density n{r). 
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Figure 5. Convergence of the energy for TV = 4, A = 1.3 and f3 = 5.0 - Panel (A) 
shows the convergence of the total energy versus the inverse number of propagators 
P~^ oc e. Shown are the results for two different choices of the parameters, a) to = 0.04, 
ai = 0.0 and b) to = 0.13, oi = 0.33, and the correct energy from CPIMC with the 
corresponding confidence interval. Panel (B) shows the decay of the average sign S 
with increasing P and panels (C) and (D) display the potential and kinetic energy V 
and K, respectively, where E = V + K. 


In Fig. 5, results are summarized for iV = 4 electrons with A = 1.3 and /d = 5, i.e., 
moderate coupling and low temperature, and panel (A) shows the convergence of the 
total energy as a function of the inverse number of propagators which is proportional to 
the imaginary time step, e oc l/P. The red diamonds [a) to = 0.04, oi = 0.0] and blue 
circles [b) to = 0.13, ai = 0.33] denote two different combinations of free parameters 
and exhibit a clearly different convergence behavior towards the exact result known from 
CPIMC, i.e., the black line. For P = 2, the energy with parameter set a) is too low by 
almost one percent. With increasing P, E increases and reaches a maximum around 
P = 6, until the curves approach the exact energy from above. For parameter set b), the 
energy converges monotonically from above and, even for P = 2, the deviation from the 
CPIMC result is as small as 0.2%. The selected energies which are listed in table 1 reveal 
that the total energy is converged for P = 15 within the statistical uncertainty. For 
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Figure 6. Influence of the relative interslice spacing to for At = 4, A = 1.3 and /3 = 5.0 
- In the left panel, the total energy is plotted versus the free parameter tg for P = 2, 
P = 3 and P = 4. The right panel shows the behavior of the average sign. 


the panels (C) and (D), the energy has been split into a potential {V) and kinetic {K) 
contribntion. For both parameter combinations, V converges monotonically, althongh 
from different directions. In addition, parameter set b) gives a mnch better resnlt 
for small P. Panel (D) reveals, that the kinetic energy K is responsible for the non- 
mo notonons convergence of E for parameter set a), which again delivers worse resnlts 
for P = 2, as compared to the bine circles. Finally, panel (B) shows the average sign S 
as a fnnction of 1/P. Both cnrves exhibit a similar decrease with an increasing nnmber 
of propagators, as it is expected. However, parameter set a) always allows for a better 
sign than b). The reason for this behavior is the free parameter to, which controls the 
relative spacing between the three time slices of an imaginary time step e. For to = 0.04, 
there are a single small and two large steps. The latter allow for more blocking, since the 
corresponding decay length in the diffnsion matrices is large as well. For to = 0.13, 
on the other hand, there are three nearly eqnal steps, each of which with a smaller decay 
length than the two large ones for parameter set a). Therefore, less blocking is possible 
and more determinants with a negative sign appear in the Markov chain. 

The different convergence behaviors of the two free parameter combinations for 
small P leads to the qnestion how to choose to and ai for optimal resnlts. To provide 
an answer, we consider the same system as in Fig. 5, and investigate the accnracy of 
the total energy as a fnnction of to, for a hxed ai = 0.33. The simnlation resnlts are 
shown in the left panel of Fig. 6 for P = 2 (red sqnares), P = 3 (bine circles) and 
P = 4 (green diamonds). All three cnrves exhibit a similar decay towards the exact 
valne starting from small to, followed by a minimnm aronnd to = 0.14 and hnally an 
increasing error for larger valnes. We note that as few as two propagators allow for an 
accnracy of |AP|/P < 2 x 10“^ for the best choice of the free parameters. Fig. 6 (B) 
shows the dependency of the average sign S on to. Again, we observe that S decreases 
with increasing to as explained dnring the discnssion of Fig. 5. In addition, it is revealed 























Permutation blocking path integral Monte Carlo 


15 




0 0.9 1.8 2.7 3.6 0 0.9 1.8 2.7 3.6 

r r 

Figure 7. Convergence of the radial density for TV = 4, A = 1.3 and /3 = 5.0 - The 
radial density n is plotted versus the distance to the center of the trap, r. In panel (A), 
the free parameters are chosen as tg = 0.13 and Oi = 0.33 and the convergence with 
P is illustrated. Panel (B) compares two different sets of free parameters, a) to = 0.13 
and Oi = 0.33 and b) to = 0.04 and ai = 0.0, for P = 2. 

that the combination of P = 4 and to = 0.01 leads to a larger sign than P = 3 and 
to > 0.10. However, the optimum free parameters allow for a higher accuracy even for 
P = 2, compared to small to with more propagators. Therefore, it turns out to be 
advantagous to use the fourth order factorization with the two free parameters despite 
the smaller average sign for the same P compared to the factorization with only a single 
daughter slice for each propagator, i.e., to = 0.0. 

Finally, we mention that the optimal choice of oi and to depends on the observable 
of interest. In Fig. 7, we investigate the effects of the free parameters on the convergence 
of the radial density distribution n{r) for the same system as in Figs. 5 and 6. The left 
panel shows u as a function of the distance to the center of the trap, r, for four different 
P and the parameter combination oi = 0.33 and to = 0.13, which has been proven to 
allow for nearly optimum energy values at P = 2, cf. Fig. 6. The black curve corresponds 
to P = 10 and is converged within statistical uncertainty. For P = 2 (red diamonds), 
there appear signihcant deviations to the latter, in particular n is too large around the 
maximum r 1.25 and too small at the boundary of the system. The P = 3 results 
(blue squares) exhibit the same trends although the differences towards the black curve 
are reduced. Finally, the density for P = 4 (green circles) cannot be distinguished from 
the converged data within the error bars. The right panel compares the density for P = 2 
with two different combinations of free parameters. The red diamonds [parameter set 
a)] correspond to the curve from the left panel and the green circles [parameter set b)] 
to oi = 0.0 and to = 0.04. The latter parameters clearly allow for a density distribution 
which is much closer to the exact results than the ai and to values which provide the 
optimal energy. 
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Figure 8. Temperature dependence for TV = 4 and A = 1.3 with to = 0-14 and 
oi = 0.33 - In the left panel, the total energy is plotted versus the inverse temperature 
/I for P = 2, P = 3 and P = 4 propagators and compared to exact CPIMC results. 
The right panel shows the behavior of the sign. 


4-2. Temperature dependence 

In the last section, we have demonstrated that the optimal choice of the free parameters 
ai and to allows for the calculation of energies with an accuracy of 0.1% with as few 
as two propagators, even at a relatively low temperature, f3 = 5.0. However, with 
decreasing T (i.e., increasing /3) the number of required propagators must be increased 
to keep the commutator error hxed. In Fig. 8, we investigate the effect of a decreasing 
temperature on the accuracy provided by a few propagators P for = 4 electrons 
at indermediate coupling, A = 1.3. The left panel shows the total energy E as a. 
function of the inverse temperature (3. We compare results for P = 2 (green circles), 
P = 3 (red diamonds) and P = 4 (blue triangles) to exact results from CPIMC (black 
stars). At larger temperature, fd < 7.0, all four datasets nearly coincide and exhibit 
the expected decrease towards the energy of the ground state. With increasing (3, the 
P = 2 results exhibit an unphysical drop because two propagators are not sufficient 
and the commutator errors become more significant. The red and blue curves exhibit 
a qualitatively similar trend, however, the energy drop is weaker and shifted to lower 
temperature. Even at 13 = 10.0, which is already very close to the ground state, three 
propagators allow for an accurate description of the system. 

In the right panel of Fig. 8, the average sign S is plotted versus the inverse 
temperature. At small (3, the wavefunctions of the electrons do not overlap and, hence, 
the system is not degenerate. With decreasing temperature, exchange effects become 
increasingly important which leads to a decrease of S. However, while for standard 
PIMC the sign is expected to exponentially decrease with (3, S seems to converge for 
PB-PIMC with P = 3 and P = 4 and exhibits an even slightly non-monotonous behavior 
for P = 2. The application of antisymmetric propagators leads to a competition with 
respect to S and (3. On the one hand, with increasing inverse temperature off-diagonal 
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matrix elements are increased, which leads to more negative determinants and, therefore, 
more negative weights in the Markov chain. On the other hand, the thermal wavelengths 

and A 2 toe are increasing with /3, which makes the blocking of large diagonal and off- 
diagonal elements more effective. Hence, the sign can even become larger with (3 once the 
system has reached the gronnd state, becanse the particle distribntion remains constant 
while more elements in the diffnsion matrix compensate each other in the determinants. 

We conclude that few propagators allow for the calculation of accurate results up 
to low temperature, {3 < 10.0. For higher /3, the system is in its ground state and finite 
temperature path integral Monte Carlo is no longer the method of choice. 

3 . 3 . Dependence on the coupling strength 

In the previous sections, we have restricted ourselves to the investigation of small systems 
to illustrate the convergence and sign behavior depending on relevant parameters. In 
this section, we demonstrate that PB-PIMC allows for the calculation of accurate results 
at parameters where no other ab initio results have been reported, so far. Fig. 9 shows 
results for N = 8 and = 20 electrons at /3 = 3.0 over a wide range of coupling 
parameters, A. In panel (A), the average sign S is plotted versus A for standard PIMC 
(squares), CPIMC (circles) and the present PB-PIMC (diamonds) with P = 2 and the 
parameter sets to = 0.14 and ai = 0.33 (A^ = 8, blue symbols) and to = 0.10 and 
Oi = 0.33 (A^ = 20, red symbols), which are known to allow for accurate energies, cf. 
Fig. 6. It is well understood that PIMC allows for the simulation of strongly coupled 
fermions, where exchange effects do not play a dominant role. With decreasing A, the 
sign exhibits a sharp drop and the sign problem prevents the simulation within feasible 
computation time for A < 2.0 and A < 5.0, respectively. Evidently, larger systems lead 
to a more severe decrease of S at larger coupling strength. CPIMC, on the other hand, 
can be interpreted as a Monte Carlo simulation on a perturbation expansion around 
the ideal quantum system, i.e., A = 0.0. Hence, the method efficiently provides exact 
results for small coupling, where the system is close to an ideal one. For A^ = 20 around 
A ~ 0.3, the sign almost instantly drops from S ~ 0.97 towards zero, and CPIMC is 
no longer applicable, without further approximation. This means that, in particular for 
larger systems, there have only been results for systems that are a) almost ideal or b) 
so strongly coupled that fermions and bosons lead to nearly equal physical properties. 
The physically particularly interesting regime where Coulomb correlations and Fermi 
statistics are significant simultaneously, has remained out of reach. 

However, the average sign from PB-PIMC exhibits a much less severe drop with 
decreasing A than standard PIMC and saturates for A < 0.7. For N = 8, the average 
sign remains above S = 0.08, which allows for good accuracy with relatively low effort. 
The small sign, S ~ 10“^, for 77 = 20 indicates that the simulations are computationally 
involved but, in contrast to PIMC and CPIMC, still feasible. In panel (B) of Fig. 9, the 
total energy E for A^ = 20 is plotted versus A over the entire coupling range and the 
statistical uncertainty from the PB-PIMC results is smaller than the size of the data 
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Figure 9. Coupling dependence for iV = 8 and iV = 20 at /3 = 3.0 - Panel (A) shows 
the average sign as a function of A for CPIMC, PIMC and PB-PIMC with N = 8 (blue 
symbols, parameter set to = 0.14 and oi = 0.33) and N = 20 (red symbols, parameter 
set to = 0.10 and ai = 0.33) and panel (B) the corresponding total energies, E, for 
the latter. In panels (C) and (D), the radial density n is plotted versus the distance 
to the center of the trap, r, for TV = 20 with A = 0.1 and A = 15.0, respectively, and 
the parameter set oi = 0.0 and to = 0.04. 


points. Both, at small and large A, the P = 2 results are in excellent agreement with 
the exact energy known from the other methods and, in addition, results are obtained 
for the particularly interesting transition region (region [II] in Fig. 1). In panel (C), 
we show the radial density for iV = 20 and low coupling, A = 0.10, calculated with 
the parameter set to = 0.04 and oi = 0.0, which has been proven effective for accurate 
densities n{r). The PB-PIMC results (red diamonds) are in excellent agreement with 
the exact CPIMC data (blue squares) over the entire r-range. For completeness, we 
mention that this combination of parameters allows for an approximately three times 
as high sign as the choice from panels (A) and (B), which was choosen to allow for a 
good energy. Panel (D) shows the density of a strongly coupled system, A = 15.0, and 
N = 20. Again, the two propagators already provide very good agreement with the 
exact curve. In Fig. 1 (B), we have shown density prohles for coupling parameters over 
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Figure 10. Influence of quantum statistics for N = 20 and /3 = 3.0 - We show the 
radial density n(r) for Fermi-, Bose- and Boltzmann-statistics in the transition region 
for A = 2.0 (A) and A = 0.7 (B). 


the entire coupling range. At A = 15 (red pluses), there are three distinct shells and 
the physical behavior is dominated by the strong Coulomb repulsion. Decreasing the 
coupling to A = 5 (green bars) leads to a reduced extension of the system, and the three 
shells exhibit a much larger overlap. At indermediate coupling, A = 2 (blue crosses), 
both the interaction and fermionic exchange govern the system. The density prohle is 
still signihcantly more extended than the ideal pendant, but n exhibits modulations 
instead of a flat curve. Decreasing the repulsion further to A = 0.7 (pink circles) leads 
to a further reduction of the extension. However, n does not approach a Gaussian- 
like prohle as for ideal boltzmannons or bosons, but continues to exhibit the density 
modulations which are characteristic for fermions. For A = 0.1, the system is almost 
ideal and the density is completely dominated by the quantum statistics. 

Finally, in Fig. 10 we compare density prohles for A^ = 20 particles at /3 = 3.0 
with Fermi-, Bose- and Boltzmann statistics. Panel (A) shows results for intermediate 
coupling, A = 2.0. The distinguishable boltzmannons (blue diamonds) exhibit a nearly 
hat prohle without any shell structure, i.e., a liquid-like behavior. The bosonic particles 
(green circles) lead to an even smoother curve, with a slightly reduced extension of the 
system. For fermions (red squares), on the other hand, the exchange already plays a 
signihcant role, as the particles exhibit an additional repulsion due to the Pauli principle, 
and n decays only at larger r. In addition, the fermionic density prohle exhibits distinct 
modulations. In panel (B), we show a comparison for smaller coupling, A = 0.7. Again, 
the boltzmannons and bosons lead to smooth density prohles which are very similar, 
despite a reduced extension of the Bose-system and an increased density around the 
center of the trap. The fermions exhibit a diherent behavior as the system is signihcantly 
more extended and the density prohle again features distinct modulations. 

In conclusion, we have presented ab initio results for the energy and the density 
for up to 20 electrons over the entire coupling range. A comparison with standard 
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Figure 11. Particle number dependence of the average sign for A = 0.1 and j3 = 3.0 
and two different combinations of the simulation parameters. 


PIMC and CPIMC has revealed excellent agreement in both the limits of weak and 
strong coupling. A more detailed investigation of the transition from the classical to the 
degenerate regime, including systematic comparisons with bosons and boltzmannons, is 
beyond the scope of this work and will be published elsewhere. 

4-4- Particle number dependence 

In the last section, we have shown that the sign problem is more severe for larger systems, 
cf. Fig. 9 (A). Here, we provide a more detailed investigation of the performance of our 
method in dependence on the particle number. In Fig. 11, the average sign S is plotted 
versus Af for A = 0.1 and fd = 3.0, i.e., a very degenerate system, with two different 
combinations of free parameters. It is revealed that S exhibits an exponential decay with 
the system size and, as usual, the smaller t^ leads to a more effective blocking. Therefore, 
the PB-PIMC approach still suffers from the fermion sign problem, and feasible system 
sizes for 2D quantum dots at weak coupling are limited to Af < 30. This is a remarkable 
result since standard PIMC simulations for A = 0.1 and (3 = 3.0 are possible only for 
Af < 4. 

5. Discussion 

In summary, we have presented a novel approach to the path integral Monte Carlo 
simulation of degenerate fermions at hnite temperature by combining a fourth-order 
factorization of the density matrix with a full antisymmetrization between all imaginary 
time slices. The latter allows to merge 3PN\ conhgurations from the standard PIMC 
formulation into a single conhguration weight, thereby efficiently grouping together 
permutations of opposite signs which leads to a signihcant relieve of the fermion sign 
problem. To efficiently run through the resulting conhguration space at arbitrary 
system parameters, we have modihed the widely used continuous space worm algorithm 
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by introducing an extended ensemble with open configurations and by temporarily 
constructing artificial trajectories. We have demonstrated the capabilities of our method 
by simulating up to = 20 electrons in a quantum dot. It has been revealed that 
the (empirical) optimal choice of the free parameters oi and to from the fourth order 
factorization allows for the calculation of energies with an accuracy of 0.1% even for 
just two propagators. For completeness, we mention that different observables lead to 
different optimal parameters. We have concluded, that it appears to be favourable to 
use two instead of a single daughter time slice for each time step e, despite the reduced 
sign for the same number of propagators. 

The investigation of the temperature dependence of the convergence with respect to 
the number of time steps P has revealed, that as few as three propagators are sufficient 
to accurately simulate fermions, up to /3 < 10.0. For larger inverse temperatures, the 
system approaches its ground state and finite temperature path integral Monte Carlo 
techniques are no longer the methods of choice. 

To demonstrate that our PB-PIMC approach allows for the calculation of accurate 
results for systems beyond the capability of any other quantum Monte Carlo technique, 
we have simulated = 20 electrons at relatively low temperature, (3 = 3.0, and 
arbitrary coupling strength. CPIMC excells at weak coupling and provides exact results 
for A < 0.3, i.e., in the region where the systems are still close to the ideal case. 
Standard PIMC, on the other hand, is applicable at strong coupling A > 5.0 where 
exchange effects are not yet dominating, until the rapid decrease of the sign renders 
any simulation unfeasible. For PB-PIMC, the sign converges for A < 0.7 and, hence, 
computations are possible at arbitrary degeneracy, in particular, in the physically most 
interesting transition region between classical and ideal quantum behavior. We find 
excellent agreement with both PIMC and CPIMC in both the limits of strong and weak 
coupling. Finally, we have demonstrated that PB-PIMC still suffers from the fermion 
sign problem, since, as expected, S exponentially decreases with the particle number. 

A possible future application of PB-PIMC to the quantum dot system might 
include the investigation of the transition from the classical to the degenerate quantum 
regime, in particular a systematic comparison of fermions to bosons and boltzmannons. 
Furthermore, it could be interesting to extend the considerations to ?)D confinements, 
e.g. [45, 46], and study the impact of quantum statistics on structural transitions [47]. 
In addition, we expect our method to be of interest for the future investigation of 
numerous Fermi systems, including the finite temperature homogeneous electron gas 
[8, 9, 10], two-component plasmas [11, 12, 13] and fermionic bilayer systems [5, 6, 7]. 
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Figure Al. Illustration of the updates Deform (left) and Swap (right) - In the 
left panel, the Deform update is executed in an open configuration. The random 
construction of an artificial trajectory (the beads marked by black arrows) is followed 
by the re-sampling of all beads between its first (start) and last (end) bead. In the 
right panel, the Swap move is demonstrated. The current head is ’connected’ to a 
random target bead on the time slice of the tail. 


Appendix A. Monte Carlo updates 


In this appendix, we present an ergodic set of Monte Carlo updates which are based on 
the usual continuous space worm algorithm [20, 40] from standard PIMC. 

(i) Deform: This update is similar to standard PIMC techniques, e.g. [43], and 
deforms a randomly constructed artificial trajectory. 

• Select a start time uniformly from all 3P slices. 

• Select a ’start’ bead on r^. 

• Select the number of beads to be changed, m G [1, M]. 

• Select m + 1 beads on the next slices according to 


7^ _ 

select 


n 

i=0 


Po(r- 


P-i-i) C 


^old 


(A.l) 


with being the normalization and the label ’old’ indicates the configuration 
before the update. 

Resample m beads in the middle according to Eq. (24): 


^ /^new ^new 
-^resample / \ ? 

Po(^0: ^tot j 

and Ctot denotes the imaginary time difference between the fixed endpoints. 
The constant M is a free parameter and can be optimized to enhance the 
performance. The update is self-balanced and the Metropolis solution for the 
acceptance probability is given by 

detpf' 


A 


Deform 


;X->X) =min l,e 


— eA# 


n 

i=0 


detp' 


,old 


(A.3) 
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with containing both the change in the potential energy and all forces. Deform 
is illustrated in the left panel of Fig. Al. 

(ii) Open/ Close: This update pair constitutes the only possibility to switch between 
open and closed configurations. The Open move is extecuted as follows: 


• Select the time slice of the new head, Thead, uniformly from all 3P slices. 

• Select the bead of the new head, rhead- 

• Select the total number of links to be erased as m G [1, M]. 

• Select m beads on the next slices from 




m— 1 

n 

i=0 


Po(l'i) Tj+i, Cj) 


the last one will be the new tail after the update. 
Delete m — 1 beads between the new head and tail. 


(A.4) 


The reverse move closes an open configuration. Let m denote the number of missing 
links between head and tail. If m > M, the update is rejected. 

• Sample m — 1 new beads according to Eq. (24) with head and tail being the 
hxed endpoints: 

imoVo(D,ri+i,ei) 




sample 


Po(^head; ^taih ^tot) 

The acceptance ratios are computed as 


m—1 
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Open 


;X ^ X) = min 1, JJ 
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Aciose(X -G X) = min | 1, 
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1 detp“' 


Ei detp° 


.old 


3CPMN 


Po(^taib ^head) ^tot) 


(A.5) 


(A.6) 


(A.7) 


The parameter /r is another degree of freedom of the algorithm and plays the same 
role as the chemical potential in the usual WA-PIMC scheme. 


(hi) Swap: The Swap move very efficiently generates exchange, i.e., allows for a switch 
between large off-diagonal or diagonal diffusion matrix elements as it is illustrated 
in the right panel of Fig. Al. Let m denote the number of missing beads between 
head and tail. 


Choose a target bead on the slice Ttaii according to 

rp _ Po(^'head) ^ti ^tot) 

target ^7; ; 

^forward 

with Sforward being the normalization. The tail itself cannot be chosen. 
Choose backwards m -I- 1 beads according to 


(A.8) 
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The head itself cannot be selected on the last slice and the last bead will be 
the new head after the update. 

’Connect’ the old head with the target bead by re-sampling the m beads 
between the slices of head and tail according to 

n 'm / new „new ^ \ 

i=oPo[^i 


^sample / \ 

POv^head: ^target? ^tot j 

The update is self-balanced and the acceptance ratio is calculated as 

detft” 


A. 


Swap 


(X -> iC) = min l,r?n 


i=0 


detp° 


old 


with the abbreviation 

—eA'H ^forward 


T] = e 


-'reverse 


(A. 10) 


(A.ll) 


(A. 12) 


and Sreverse being the normalization of the selection of the target bead from the 
reverse move. 


(iv) Advance/ Recede: These updates move the head forward (backward) in the 
imaginary time. However, they are optional and, in principle, not needed for 
ergodicity. The Advance move is executed as follows: 


Calculate the number of missing beads between head and tail, a. If a = 0, the 
update is rejected. 

Select the number of new beads to be sampled, m G [1, a]. 

Sample the position of the new head from po(i'head, ^tot)- 

Sample the m — 1 beads between old and new head according to Eq. (24) 


nr=^po(rr 


- sample 


^i+l 


Pi) 


Po(rhead,rSietot) 


(A.13) 


The reverse move is given by Recede. Let k denote the total number of beads which 
can be removed. If k = 0, the update is rejected. 


• Select the total number of beads to be removed as m G [1, k]. 

• Select m beads backwards starting from the old head from 


Eselect 


m—1 

n 


Po(rr 


^ j+1)' 


(A.14) 


1=0 ‘ 

with being the normalization. The last one will be the new head after 
the update. Here “new” denotes new with respect to Advance, since the 
coordinates are pre-existing for the Recede move. Delete the m beads between 
the new head and tail. 


This gives the acceptance ratios 

^Advance (X -A X) = miu 
^Recede (X ^ X) = miu 



(A.15) 
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with the definition 

e = . (A. 16 ) 

K 

The presented list of Monte Carlo moves constitntes an ergodic set of local npdates, 
which allows for an efficient sampling of both the extended confignration space and a 
canonical Markov chain. 
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